% -------------------------------------------------------------------------
%                   Price Discrimination in Selection Markets
%                                   Andre Veiga
%                   Review of Economics and Statistics, 2023
% -------------------------------------------------------------------------




%% setup

rng('default'); % set the seed each time
clear;
beep on;
warning('off','all');

% path where stuff is saved
output_path = 'output/';

% set(groot, 'defaultLegendInterpreter','latex');
set(groot,'defaultTextInterpreter','latex');


%% make one graph illustrating equilibrium

% make one graph as an example of how the curves look and what
% equilibrium is
option = 0;
[uA, cA, uB, cB ] = gen_u_c(option);
[rows,cols,pp] = deal(1,1,0); clf
pp=pp+1; subplot(rows,cols,pp)
make_mkt_graph(uA, cA,  'Equilibrium', 0);
mysave( gcf, output_path, 'eql_eg', 1*[rows,cols]);


%% make graphs illustrating price discrimination in various settings

tic_all = tic;
MODELS = [ 1:7 ]; % there are 7 sets of markets to consider
for i= 1:numel(MODELS)
    option = MODELS(i);
    [uA, cA, uB, cB ] = gen_u_c(option); % generate cost and WTP for both markets
    file_name_i = strcat('small_', num2str(option) );
    make_delta_graphs(uA, cA, uB, cB, file_name_i, output_path);
end
toc_all = toc(tic_all);
fprintf('Everything ran in %4.2fs \n', toc_all)






%% auxiliary functions










%% check types

function [z, p_z] = compute_z(u, c)

% compute z = -Q'/Q;

assert( size(u, 1)==1 , 'u must be row vec')
assert( size(c, 1)==1, 'c must be row vec' )

% evaluate at a fixed number of points
nr_points_to_eval = 20;
idx = round( linspace(1, numel(u), nr_points_to_eval));
p_z = u(idx);
cc = c(idx);

q = NaN(size(p_z));
% compute quantity at each price
for i=1:numel(p_z)
    q(i) = mean( u >= p_z(i) );
end
% compute z
deriv_q = gradient(q)./gradient(p_z);
sigma = - deriv_q./q;
assert( isequal(size(sigma), size(q)), 'sigma and q should have same size')
z = sigma .* (p_z - cc);

% checks on outputs
assert( size(z,1)==1 )

end



function [] = check_u_c( u, c )

% checks that types u,c are OK

% both markets should have the same number of types (same market size)
% and WTP uA, uB should be increasing
% all of these should be row vecs
assert( isequal(numel(u), numel(c)), 'u,c must have same number of elements' )
assert( all(diff(u)>=0), 'u must be increasing')
assert( size(u, 1)==1 , 'u must be row vec')
assert( size(c, 1)==1, 'c must be row vec' )
assert( all(c>=0) , 'must have c>=0')
assert( all(u>=0) , 'must have u>=0')


% check that z is increasing
z = compute_z(u,c);
condition_1 = all( diff(z) >=0 );
if ~condition_1
    disp('z(p) not globally increasing!')
end

% check that c crosses u once and from above
% check most eager type has (weakly )positive surplus
surplus = u-c;
assert( surplus(end)>=0, 'most eager type must have u-c>=0')


% check that surplus is always positive or crosses 45 degree line once and
% from above
% check MC cuts the 45 degree line only once and from above, the surplus
% goes from negative (at low WTP) to positive (at high WTP), so the sign
% change is positive
% it's important that this accounts for equality as well, since it's OK to have u=c
if all(surplus >= 0) 
    return
else
    assert( surplus(1)<0 )
    diff_sign_surplus = diff(sign(surplus));
    sign_changes = diff_sign_surplus( diff_sign_surplus~=0);
    assert( numel(sign_changes)==1, 'c(p) can only cross u once' ) % if it does not strictly cross, this will have been caught by the previous if condition
    assert( sign_changes>=0 , 'surplus must go from negative to positive')
end


end





%% generate types

function [uA, cA, uB, cB ] = gen_u_c(option)

% generates WTP and costs for a given market under various scenarios,
% determined by the input option
% each option includes a brief description of the characteristics of that
% pair of markets in the variable descript

fprintf('\nGenerating Individuals for option %u... \n', option)

% this is the default number of simulated individuals but in some instances
% there can be more or fewer indivs simulated
nr_indivs = 2000;

switch option
    
    case 0
        descript = 'illustration of equilibrium graph';
        uA = linspace(0, 200, nr_indivs);
        cA = uA/2 + 20 ;
        [uB, cB] = deal(uA, cA); % this is just so the checks at the end of this function can run smoothly.
        
    case 1
        descript = 'optimal policy is interior';
        uB = linspace(0, 100, nr_indivs);
        cB = uB/4 ;
        uA = linspace(0, 200, nr_indivs);
        cA = uA/2 ;
        
    case 2
        descript = 'optimal policy is full PD';
        uB = linspace(0, 100, nr_indivs);
        cB = uB/2 + 10 ;
        uA = linspace(0, 100, nr_indivs);
        cA = uA/3 + 50  ;
        
    case 3
        % here it's important that there are many simulated individuals for the
        % graph to be smooth
        descript = 'optimal policy is full CR';
        nr_indivs = 10000;
        uB = linspace(0, 50, nr_indivs);
        cB = uB*0.02 + 2 ;
        uA = linspace(0, 50, nr_indivs);
        cA = uA*(0.1) + 1 ;
        
    case 4
        % for some delta, pB is the monopoly
        % price in market B; the optimal delta is higher than the one at
        % which market B collapses
        descript = 'market B collapses at full CR';
        uB = linspace(0, 100, nr_indivs);
        cB = uB/3 +1 ;
        uA = linspace(0, 200, nr_indivs) + 60;
        factor = 1.2;
        cA = uA*factor - uA(end)*(factor-1) - 20 ;
        
        
    case 5
        descript = 'market A collapses at full PD, interior policy is optimal';
        uB = linspace(0, 200, nr_indivs) + 50;
        cB = uB/3 + 50 ;
        uA = linspace(0, 200, nr_indivs) + 60;
        factor = 1.25; % the bigger is the factor, the more MC is below WTP
        cA = uA*factor - uA(end)*(factor-1);
        
        
    case 6
        descript = 'market A collapses at full PD, full PD is optimal';
        uB = linspace(0, 200, nr_indivs) + 50;
        cB = uB/2 + 50 ;
        uA = linspace(0, 200, nr_indivs) + 60;
        factor = 1.25; % the bigger is the factor, the more MC is below WTP
        cA = uA*factor - uA(end)*(factor-1);
        
        
        
    case 7
        % the HHW market
        descript = 'HHW markets';
        
        % the levels of generosity. I assume the outside option gives zero
        % insurance, so xL = 0
        xH = 0.8;
        xL = 0;
        
        % sqrt_nr_indivs is square root of number of individuals to generate (types are
        % generated using quadrature in 2D, so sqrt_nr_indivs^2 individuals
        % end up being generated in each market
        sqrt_nr_indivs = 150;
        
        % generates mu (expected risk) and v (insurance value) based on
        % estimates from HHW, for each of 8 age-based markets
        [ MU, VV ] = simm_consumers_HHW( sqrt_nr_indivs );
        
        markets_A = 5:8; % the high cost markets, the oldest
        markets_B = 1:4; % the low cost markets, the youngest
        
        % to combine various age-based, stack the vectors of types
        mu_A_HHW    = make_col_vec( MU(:, markets_A)  );
        v_A_HHW     = make_col_vec( VV(:, markets_A)  );
        mu_B    = make_col_vec( MU(:, markets_B)  );
        v_B     = make_col_vec( VV(:, markets_B)  );
        
        % compute cost and WTP for market A, based on mu and v
        % since xL=0, cost is equivalent to x*mu
        cA = (xH-xL)*mu_A_HHW;
        uA = cA + ( gg(xH) - gg(xL) )*v_A_HHW;
        
        % compute cost and WTP for market B
        cB = (xH-xL)*mu_B;
        uB = cB+ ( gg(xH) - gg(xL) )*v_B;
        
        % sort types and cost so that WTP is increasing
        [uA, idx] = sort(uA);
        cA = cA(idx);
        [uB, idx] = sort(uB);
        cB = cB(idx);
        
        % outputs should be row vecs
        [uA, cA, uB, cB] = deal(uA', cA', uB', cB');
        
        
    otherwise
        error('invalid value for option')
        
end


% checks on outputs
check_u_c( uA, cA)
check_u_c( uB, cB )

% output message
fprintf('Generated %u Individuals (per market) for option %u (%s) \n', numel(uA), option, descript)


end



%% make delta graphs


function [] = make_delta_graphs(uA, cA, uB, cB, file_name, save_folder)

% makes the graphs of prices and total welfare as delta (the constraint on price
% differences) changes


% find eql prices under full PD for each market
[ p_star_A ] = find_eql_price(uA, cA);
[ p_star_B ] = find_eql_price(uB, cB);

% build grid of deltas
delta_max = p_star_A - p_star_B;
nr_deltas = 30;            % number of points in the grid of deltas
DELTAS = linspace(0, delta_max, nr_deltas);

% pre-allocate memory for outputs
[ qA, qB, pA, pB, piA, piB, wA, wB ] = deal(NaN(nr_deltas,1));

% compute outcomes for each delta
for i = nr_deltas: - 1: 1
    delta_i = DELTAS(i);
    
    % find equilibrium prices given the types and the level of the constraint delta_i
    [ pa_i, pb_i ] = find_eql_price_constrained(uA, cA, uB, cB, delta_i);
    
    % compute outcomes in each market for the equilibrium prices: quantity,
    % profit and welfare
    [ qA(i), ~, piA(i), wA(i) ] = outcomes(uA, cA, pa_i);
    [ qB(i), ~, piB(i), wB(i) ] = outcomes(uB, cB, pb_i);
    [ pA(i), pB(i) ] = deal(pa_i, pb_i);
end


% make graphs
[rows,cols,pp] = deal(1,4,0); clf 

% set graph_q=1 to make the x-axis quantity (otherwise x-axis is price)
graph_q = 0;

% graphs WTP, AC, c for markets A and B
pp=pp+1; subplot(rows,cols,pp)
[ pA_monop ] = make_mkt_graph(uA, cA, 'Market A', graph_q);

pp=pp+1; subplot(rows,cols,pp)
[ pB_monop ] = make_mkt_graph(uB, cB, 'Market B', graph_q);


% plot equilibrium price in each market as a function of the constraint delta
pp=pp+1; subplot(rows,cols,pp)
plot(DELTAS, pA, '-o'); hold on
plot(DELTAS, pB, '-s');
yline(pA_monop, ':b' )
yline(pB_monop, '--r')
legend('$p_A$', '$p_B$', '$\hat{p_A}$', '$\hat{p_B}$',  'Location', 'Best', 'AutoUpdate','off', 'Interpreter','latex')
xlabel('$\delta$')
title(['$p_A, p_B, \Delta=$' num2str(delta_max, '%.2f')])


% graph total welfare as a function of the constraint delta
pp=pp+1; subplot(rows,cols,pp)
w_tot = wA + wB;
plot( DELTAS, w_tot, '-o' ); hold on
[max_w, idx] = max(w_tot);
delta_star = DELTAS(idx);
scatter(delta_star, max_w, 150, 'd', 'filled')
xlabel('$\delta$')
title(['Total Welfare, $\delta^*=$' num2str(delta_star, '%.2f')])

% save graph
mysave( gcf, save_folder, file_name, 1*[rows,cols]);


% what is the welfare gain relative to full PD?
w_PD = w_tot(end);
gain_in_welf = max(w_tot) - w_PD;
gain_in_welf_pct = 100*(max(w_tot) - w_PD)/abs(w_PD);
fprintf('welfare gain (at optimum) relative to full PD: %4.2f $ = %4.2f %% \n', gain_in_welf, gain_in_welf_pct)
fprintf('delta opt= %4.2f, pA*-pB* = %4.2f, ratio = %4.2f \n', delta_star, delta_max, delta_star/delta_max)

end




%% risk aversion

function [y] = gg(x)
% this is the expression for risk adjustment
% utility is u = x*mu + v*gg(x)-p
y = (1/2)*( 1-(1-x).^2);
end







%% make graph for a single market market

function [ pA_monop ] = make_mkt_graph(uu, cc, mkt_label, graph_q)

% creates the graph for a single market, showing WTP, AC and c
% if graph_q=1, the x-axis is quantity
% mkt_label is the title of the graph

% starting message
fprintf('Making market graph for %s \n', mkt_label )


nr_points = numel(uu);
ppp = [1:nr_points]'; % default: use every value in uA, cA;
max_nr_points_p = 15;
if nr_points>=max_nr_points_p
    % if there are too many values in uA, cA pick a uniformely distributed subset of uA, cA, etc at which to
    % evalute all functions
    ppp = round( linspace(1, nr_points, max_nr_points_p) )';
end

% find equilibrium price for market A
[ p_star_A ] = find_eql_price(uu, cc);

% find monopoly price for market A
ppA = uu; % the grid of prices to try = the grid of WTP
[ q, AC, pi, ~, MC ] = outcomes(uu, cc, ppA);
[~, idx ] = max(pi);
pA_monop = ppA(idx);


switch graph_q
    
    case 0
        % x-axis is prices; evaluated at the points ppp (which are a
        % uniform subset of uA, cA
        plot( ppA(ppp), ppA(ppp), '-o'); hold on
        plot( ppA(ppp), AC(ppp), '-s')
        plot( uu(ppp), MC(ppp), '-d' )
        xlabel('p ')
        
    case 1
        % x-axis is quantity
        plot( q(ppp), ppA(ppp), '-o'); hold on
        plot( q(ppp), AC(ppp), '-s')
        plot( q(ppp), MC(ppp), '-d' )
        xlabel('q')
        
    otherwise
        error('invalid value for graph_q')
end

legend('p', 'AC', 'c', 'Location', 'Best', 'AutoUpdate', 'off', 'Interpreter','latex')
title([mkt_label, '$, p^*=$' num2str(p_star_A, '%.2f')])


end




%% solve for prices under the constraint


function [y] = sum_profits( uA, cA, uB, cB, delta, pb )

% computes the sum of profits across the two markets for a
% given value of the price in market B, pb, and a value of the constraint
% delta, so that pa=pb + delta
% at the equilibrium value of pB, this function must vanish

[ ~, ~, piB ] = outcomes(uB, cB, pb);
pa = pb + delta;
[ ~, ~, piA ] = outcomes(uA, cA, pa);
y = piA + piB;

end




%% find equilibrium prices when prices are constrained by regulation



function [ pA, pB ] = find_eql_price_constrained(uA, cA, uB, cB, delta)

% this function finds the equilibrium price pB, for a given value of the
% constraint delta
% 1) computes the full PD prices for each market
% 2) if the delta > difference in full PD prices, then equilibrium is full PD
% 3) otherwise find the value of pb that makes the sum of profits vanish
% If sum_pi <=0 for all values of pb, then there is market collapse


% find eql prices in each market
[ p_star_A ] = find_eql_price(uA, cA);
[ p_star_B ] = find_eql_price(uB, cB);

% check that market A is the strong market
assert( p_star_A > p_star_B, 'market A should be the strong market')

% if delta is too large, then we are at full PD, so equilibrium has been
% found
if abs(p_star_A - p_star_B)<=delta
    pA = p_star_A;
    pB = p_star_B;
    return
end



% function that finds equilibrium; at equilibrium, the sum of profits must
% vanish
sum_profit_pb = @(pb) sum_profits( uA, cA, uB, cB, delta, pb );



% number of prices pB to check for sum_pi=0; this is something that slows
% the code down somewhat
grid_p_size = 1e4;
PPP = linspace(0, max([uB, uA]), grid_p_size);
eval_sum_profit = sum_profit_pb(PPP);

% check if the sum of profit crosses zero (ie, if it ever changes sign)
% if mkt_B_collapse==1, then market B collapses
mkt_B_collapse = all( diff( sign(eval_sum_profit)) ==0 );

switch mkt_B_collapse
    case 1
        fprintf('find_eql_price_constrained, delta= %4.2f : no pB st (sum pi)=0, setting pA=pA*,  pb=pA*-delta. \n', delta)
        pB = p_star_A - delta;
        
    case 0
        idx = find(diff(eval_sum_profit>=0),1); % this is the price where sum_pi changes sign for the FIRST time
        pB_starting_val = PPP(idx);
        
        % fine tune the equilibrium value of price using fzero; this usually doesn't make a big difference but
        % it does for some markets where the max of delta is small
        options = optimset('Display', 'off' );
        [ pB, ~ ] = fzero(sum_profit_pb, pB_starting_val, options);
        
        % if the result is NaN, try a starting value at the very top of the
        % WTP distribution
        if isnan(pB)
            [ pB, ~ ] = fzero(sum_profit_pb, max(uB), options);
        end
        
        % if the result is still NaN, then allocate the starting value;
        % this should never happen and if it does a warning is printed to
        % the command window
        if isnan(pB)
            fprintf('find_eql_price_constrained: (pb after fzero)=NaN, using pB=pB_starting_val! \n')
            pB = pB_starting_val;
        end
        
    otherwise
        error('invalid value of mkt_B_collapse')
end

sum_profit_at_eql = sum_profit_pb(pB);
pA = pB + delta;

do_graph = 0;
% set do_graph = 1 to produce a graph that illustrates where the
% equilibrium is; this slows down the code somewhat
if do_graph == 1
    
    pause_lenght = 0;
    [rows,cols,pp] = deal(1,1,0); clf
    
    pp=pp+1; subplot(rows,cols,pp); hold on
    plot(PPP, eval_sum_profit);
    yline(0)
    scatter(pB, sum_profit_at_eql, 100, 'd', 'filled')
    legend('sum \pi', '0', 'eql',   'Location', 'Best','AutoUpdate','off')
    xlabel('p')
    xline(pB)
    
    if mkt_B_collapse==1
        title(['no p_B st (sum \pi)=0, setting p_B=p_A*-\delta'])
        pause(10)
    else
        title([ '\delta=' num2str(delta, '%4.2f') ', p_B*=' num2str(pB, '%4.2f')  ])
    end
    
    pause( pause_lenght )
end

end





%% computes the margin p-AC

function [y] = compute_margin(u, c, p)

% computes the margin p-AC, for a given market at a given vector of prices

% checks on inputs
% p = make_col_vec(p);
assert( size(p,2)==1 )

% compute AC and margin
[ ~, AC, ~ ] = outcomes(u, c, p);
y = p-AC;

% checks on outputs
assert( size(y,2)==1 )

end



%% find equilibrium price


function [ eql_price ] = find_eql_price(u1, c1)

% this function computes the equilibrium price for a single market (ie, the full PD price for that market)
% this is the lowest price at which p=AC


% compute the margin on a grid of prices, this yields a starting value for
% the search and also checks if the margin is everywhere positive or
% everywhere negative
% if there are multiple points where p=AC, min() selects the
% first one which will be the lowest price, and that is indeed the
% equilibrium
size_p_grid = 1e3;
PP = linspace( min(u1), max(u1), size_p_grid)';
margin = compute_margin(u1,c1, PP);
[~, idx] = min(abs(margin));
x0 = PP(idx);

% check for extreme cases where q=1 or q=0
if all( margin>0 )
    eql_price = 0;
    disp(['find_eql_price: p>AC(p) for all p, p*=0'])
    return
    
elseif all( margin<0 )
    eql_price = max(u1);
    disp(['find_eql_price: market collapse, p<AC(p) for all p, set p*=max(u)'])
    return
end



% fine-tune the above value of equilibrium price found from the grid, using
% fzero
% use x0 from above as the starting value, and fine-tune it using fzero
fn_find_zero = @(p) compute_margin( u1, c1, p );
options = optimset('Display','off' );
eql_price = fzero(fn_find_zero, x0, options);


end



%% outcomes

function [ q, AC, pi, CS, MC ] = outcomes(u1, c1, PP)

% this function computes outcomes for a given distribution of types u,c and
% a given vector of prices PP

if isnan(PP)
    [ q, AC, pi, CS ] = deal(NaN);
    fprintf('outcomes: PP=NaN \n')
    return
end


% make sure inputs are column vectors
PP = make_col_vec(PP);
u1 = make_col_vec(u1);
c1 = make_col_vec(c1);
nr_points = numel(PP);

% pre-allocate memory for outputs
[ q, AC, pi, CS ] = deal( NaN(nr_points, 1));

method = 1;
switch method
    
    case 1
        % turns out this method with the loop is faster, probably because
        % it involves a number of means conditional on subsets, and
        % considering these smaller vectors speeds things up
        for i = 1:nr_points
            p_i = PP(i);
            surplus = u1 - p_i;
            buying_probs = u1 >= p_i;
            q(i) = mean( buying_probs );
            AC(i) = mean( c1( buying_probs ) );
            CS(i) = q(i) * mean( surplus (buying_probs)  );
            pi(i)  = q(i) *( p_i - AC(i) );
            
            % if nobody buys, then AC is the MC of the most eager
            % individual, and CS = pi =0
            if q(i)==0; [ AC(i), CS(i), pi(i) ] = deal( c1(end), 0, 0); end
        end
        
    case 2
        % rows are values of u,c; columns are values of P
        PP = PP';
        buy = u1 > PP;
        
        q = mean( buy, 1);
        CC = repmat(c1, 1, nr_points);
        assert( isequal(size(CC), size(buy)) )
        AC = mean( CC.*buy, 1) ./ q;
        CS = mean( (u1 - PP).*buy, 1) ./ q;
        pi = q .* (PP - AC);
        
        % if nobody buys, set these to zero
        AC( q==0 ) = c1(end);
        CS( q==0 ) = 0;
        pi( q==0 ) = 0;
        
        % turn everything into col vecs
        [q, AC, CS, pi] = deal(q', AC', CS', pi');
        
        assert( isequal(size(q), [nr_points,1]) , 'q should be a column vector')
        assert( isequal(size(AC), [nr_points,1]) , 'AC should be a column vector')
        assert( isequal(size(pi), [nr_points,1]) , 'pi should be a column vector')
        assert( isequal(size(CS), [nr_points,1]) , 'CS should be a column vector')
        
    otherwise
        error('no such method available')
end

% compute total cost and then use this to compute marginal cost MC
tot_cost = q.*AC;             % total cost
MC = gradient(tot_cost)./gradient(q);

% check on outputs
assert( all(CS>=0), 'CS cannot be negative')


end



%% generate consumers according to the moments estimated by Handel, Hendel and Whinston ECMA 2015 (HHW)

function [ MU, VV ] = simm_consumers_HHW( sqrt_n_indivs )

% simulates consumer according to the HHW estimates

mean_v = 6.7e4;
var_v = 9.8e9;
rho = 0.56;

% distrubtion of health expenditures from HHW paper, Table III (page 1285)
MEAN_mu     = [ 3112, 3766, 4219, 5076, 6370, 7394, 9175, 10236];
VAR_mu      = [ 4918, 5473, 5304, 5942, 6874, 7116, 7414, 7619].^2;
nr_markets = numel(MEAN_mu);


% chosing the mean of v
MEAN_v = mean_v*ones(1,nr_markets);
VAR_v = var_v*ones(1,nr_markets);
RHO = rho*ones(1,nr_markets);

% moments of the distribution of WTP and cost
MOMENTS = [MEAN_mu' , VAR_mu' , MEAN_v', VAR_v' , RHO' ];

% pre-allocate memory
[MU, VV] = deal(NaN( sqrt_n_indivs^2, nr_markets));

for i=1:nr_markets
    % the function biv_logn takes as inputs the means, variances, etc of
    % the underlying gaussian variables; the function
    % recover_gauss_from_lognormal computes these gaussian moments based on
    % the desired lognormal moments
    lognormal_moments_i = MOMENTS(i,:);
    [ gaussian_moments_i ] = recover_gauss_from_lognormal( lognormal_moments_i );
    [ mean_log_mu_i, var_log_mu_i, mean_log_v_i, var_log_v_i, log_rho_i ] = deal_cols( gaussian_moments_i );
    
    [mu_i, v_i] = biv_logn( mean_log_mu_i, var_log_mu_i, mean_log_v_i, var_log_v_i, log_rho_i, sqrt_n_indivs);
    MU(:,i) = mu_i;
    VV(:,i) = v_i;
end


end


